前排提醒：如果数学公式看不清楚，可以把鼠标悬停上去，会放大。
简要题意
给定一张 $n$ 个点的简单无向图，每个点可能有一个颜色 $c_i > 0$, 也可能没有颜色，记作 $c_i = 0$.
你需要找出尽可能多的简单路径，使得：
- 这些路径间两两不交（不经过重复点）；
- 路径的起点与终点需要是有颜色但是颜色不同的点；
- 除了起点终点外，路径只能经过没有颜色的点。
需要构造出答案。
$n \leqslant 300$.
背景与参考文献
这个问题是经典问题，叫做 S-path problem。
本文将陈述此问题的 $O(n^3)$ 算法。若使用更快的 $O(n^w)$ 的矩阵乘法，则可以优化到 $O(n^w)$。若将来 $w=2$，则复杂度为 $O(n^2\log n)$。
参考文献为
理论与算法
1. 问题转化
我们先把原问题转化为：
问题一：选出原图中尽量多的边，使得这些边构成森林，且森林中每棵树里最多有两个有颜色的点，并且不能有两个相同颜色的点。
这样的话，最优解中每棵树里肯定要么有一个有颜色的点，要么有两个有颜色且颜色不同的点（除非原图中某个连通块里一个有颜色的点都没有，这样的连通块显然可以直接扔掉）。
所以原问题的答案一定就是总的有颜色的点的个数，减去新问题中连通块个数。而新问题中连通块个数又等于 $n$ 减去选出的边数。
因此选的边越多，得出的简单路径就越多。这些简单路径也可以由新问题里选出的边轻易构造出来。
这一步转化是容易的，我们接下来会将这个新问题再转化为一个线性代数问题。
2. 线性拟阵配对问题
这个线性代数问题如下。
问题二(线性拟阵配对)：
给定 $N$ 维空间里的 $M$ 对向量，即 $2M$ 个向量 $(b_1, c_1), (b_2, c_2), \dots, (b_M, c_M)$.
你需要找出其中尽可能多对，使得这些对向量线性无关。
如果你并不知道“线性无关”的定义，建议参考OI-Wiki: 线性基。
此外，OI-Wiki 中线性代数相关的其他概念也可能对你理解本文有帮助。
如何将前面的问题一转化为该问题呢？有一个非常精妙的构造如下。
我们先取一个大空间 $\mathbb{R}^{2n}$，其中每个向量 $x$ 的坐标记作 $x_1, x_2, \dots, x_{2n}$。
此外我们记 $e_i$ 为第 $i$ 个分量是 $1$, 别的分量是 $0$ 的向量。
为了方便，我们记 $x(u) = (x_{2u-1}, x_{2u})$ 为一个二维向量。
接下来对原图中每条边 $(u, v)$，我们将其对应一个二维子空间
\begin{equation}
L_{u,v} = \{ x \in \mathbb{R}^{2n} \mid x(u) + x(v) = 0 \text{ 且 } \forall w \notin \{u,v\}, x(w) = 0 \}
\end{equation}
也可以说是对应了两个向量 $b_{u,v} = e_{2u-1} - e_{2v-1}, c_{u,v} = e_{2u} - e_{2v}$。
如果只是这样的话，那么只有当选出的边构成一个环的时候，才会线性相关：
比如如果选出了 $x - y - z - x$ 这样三条边，那么 $b_{x,y} + b_{y,z} + b_{z,x} = 0, c_{x,y} + c_{y,z} + c_{z,x} = 0$。
为了满足我们的额外条件（每棵树最多有两个有颜色的点，并且颜色不能相同），还需要引入一个额外的东西。
对于每个颜色，我们选取 $\mathbb{R}^2$ 里的一个非零向量 $f_c$，并且他们两两不平行（即不存在 $c \neq c'$ 使得 $f_c$ 是 $f_{c'}$ 的若干倍）。
此外记 $f_0 = 0$。
对每个点 $u$，我们取一个向量 $a_u$，使得 $a_u(u) = f_{c_u}$，而 $a_u(v) = 0$.
也就是说 $a_u$ 的 $u$ 分量是二维向量 $f_u$，其他分量是 $0$。如果 $c_u = 0$ 那么也有 $a_u = 0$.
可以发现，如果我们选出了一些边，这些边连通了某两个点 $s, t$，那么就可以用这些边对应的向量加出 $L_{s, t}$ 的元素，即使 $(s, t)$ 这条边不存在或者没被选择。
那么如果我们选出的边连接了两个颜色相同的点，比如 $1$ 和 $2$，就有 $a_1 - a_2 \in L_{1, 2}$。
而如果我们选出的边连接了三个有颜色的点，比如 $1, 2, 3$ 且其颜色也是 $1, 2, 3$，
我们可以找出三个实数 $x_1, x_2, x_3$ 使得 $x_1 f_1 + x_2 f_2 + x_3 f_3 = 0$。那么总可以找出 $L_{1, 2}$ 和 $L_{1, 3}$ 里的元素，加起来恰好得到
$x_1 f_1 + x_2 f_2 + x_3 f_3$。
也就是说，如果我们强制选择 $f_1, \dots, f_n$ 这些向量（除了其中的零向量），接下来再选择一些边（以及他们对应的若干对向量 $(b_{u,v},c_{u,v})$），
那么选出的所有向量线性无关当且仅当这些边满足问题一中的条件。
怎么强制选择这些向量呢？可以考虑把所有 $b, c$ 向量对都投影到某个和他们正交的子空间上，这样就相当于已经选择这些向量了。
3. 转化的具体实现
首先，我们先（假装）求出这 $2m$ 个向量 $b_{u,v} = e_{2u-1} - e_{2v-1}, c_{u,v} = e_{2u} - e_{2v-1}$.
接下来我们做投影。注意我们其实不需要正交，只需要把他们投影到某个子空间 $W$ 上，满足 $\mathbb{R}^{2n} = W \oplus \mathrm{Span}(f_1, \dots, f_n)$（直和）即可。如果用 OIer 比较熟悉的说法，就是用这些 $f$ 对他们做高斯消元。
所以我们不妨设 $x_c = (1, -c)$，则 $f_u = e_{2u - 1} - c_u e_{2u}$。
用 $f_u$ 对某个向量做高斯消元的话，就是把这个向量的第 $2u$ 项系数加上 $c_u$ 倍的 $2u - 1$ 项系数，然后把第 $2u - 1$ 项系数置为 $0$。
也就是
v[2*u] += c[u] * v[2*u - 1]; v[2*u-1] = 0;
（然后我们也可以把第 $2u-1$ 维直接删掉，反正所有向量的这一维都是 $0$）
4. 新问题的求解
接下来我们就要求解这个线性代数问题。
我们先定义一个记号：若 $S, T \subseteq \{ 1, \dots, n \}$ 是两个子集，$M$ 是 $n \times n$ 的矩阵，那么 $M_{S, T}$ 是只取行序号属于 $S$ 列序号属于 $T$ 的位置构成的 $|S| \times |T|$ 子矩阵。如果 $S = $ 全集或者 $T = $ 全集，我们也简记作 $M_{*,T}, M_{S,*}$.
我们给出本文中四个不加证明的定理：
首先是解决线性拟阵配对问题的核心定理：
定理 [Lovász]：
给定 $N$ 维空间里的 $M$ 对列向量，即 $2M$ 个向量 $(b_1, c_1), (b_2, c_2), \dots, (b_M, c_M)$.
我们引入 $M$ 个未知数 $x_1, \dots, x_M$，并且在他们的分式域里考虑矩阵（即矩阵的每个分量可以是这些未知数的有理式）。
定义矩阵
\begin{equation} Y = \sum_{i = 1}^M x_i (b_i \wedge c_i). \end{equation}
其中 $b_i \wedge c_i = b_i c_i^T - c_i b_i^T$。
则 $\mathrm{rank}(Y) = 2k$，其中 $k$ 是问题二的答案，即最多可以选出的线性无关向量对数。
$\wedge$ 读作 wedge。
其次是让我们摆脱掉矩阵分量里的未知数的引理：
定理 [Schwartz-Zippel 引理]：
在模素数 $p$ 的意义下，若 $f(x_1, \dots, x_n)$ 是不超过 $d$ 次的非零多元多项式，那么任取 $x_1, \dots, x_n$ 为 $0, \dots, p-1$ 的随机值，
$f(x_1, \dots, x_n) = 0$ 的概率不超过 $\frac{d}{p}$.
然后是让我们之后算法的重点：
定理 [Harvey]：
设 $M$ 是一个 $n \times n$ 可逆矩阵, $N = M^{-1}$。若 $\tilde{M}$ 和 $M$ 长得很像——具体地说，有一个子集 $S$，
使得除了 $\tilde{M}_{S,S} \neq M_{S,S}$ 之外 $\tilde{M}$ 的分量都和 $M$ 相等。记 $\Delta = \tilde{M} - M$。那么
- $\tilde{M}$ 可逆当且仅当 $I + \Delta_{S,S}N_{S,S}$ 可逆。
- 若第一条成立，则 $\tilde{M}^{-1} = N - N_{*,S}(I + \Delta_{S,S}N_{S,S})^{-1}\Delta_{S,S}N_{S,*}$;
- 若第一条成立，则 $(\tilde{M}^{-1})_{T,T} = N_{T,T} - N_{T,S}(I + \Delta_{S,S}N_{S,S})^{-1}\Delta_{S,S}N_{S,T}$。
最后是其实可能不必要的一个玩意。
定理：
若 $A$ 是斜对称矩阵，即 $A^T = -A$，而 $\mathrm{rank}(A) = r$，则存在 $|S| = r$ 使得 $A_{S, S}$ 可逆。
具体地，只要 $A$ 中 $S$ 对应的行线性无关（即 $\mathrm{rank}(A_{S,*}) = r$），就有 $A_{S,S}$ 可逆。
如果没有斜对称的条件，那么只能知道存在 $|S| = |T| = r$ 使得 $A_{S,T}$可逆。
对可能不熟悉的读者：秩 $\mathrm{rank}$ 就是 Gauss 消元后剩下的非零行数。$n \times n$ 的矩阵 $A$ 可逆等价于 $\mathrm{rank} A = n$。
首先，我们可以取一个足够大的素数 $p$，比如 $998244353$，做模数。然后我们随机选 $m$ 个数 $x_1, \dots, x_m$ 当作“未知数”。
根据 Schwartz-Zippel 引理，只要我们接下来处理的都是多项式（比如，算矩阵的秩相当于让某个子矩阵的行列式非 $0$），就基本上不会算错。
接下来，我们算出前面这个矩阵 $Y = \sum_{i = 1}^M x_i (b_i \wedge c_i)$，并计算它的秩。但是光算出这个还不行，我们还需要找到方案。
为了方便，我们可以求出 $Y$ 的最大行向量线性无关组 $S_0$，
那么由于 $Y$ 斜对称，就知道 $Y_{S_0, S_0}$ 可逆，因此只需要关心这些行列即可。
也就是说接下来我们把所有 $b_i, c_i$ 向量里其他分量都扔掉，只留下 $S_0$ 对应的这些维分量。
这样的话 $Y$ 就可逆，我们可以算出他的逆 $N = Y^{-1}$。
我们找到方案的办法非常简单粗暴：对于每一对 $(b_i, c_i)$，我们把 $x_i (b_i \wedge c_i)$ 从 $Y$ 里减去，看看它是否还可逆。如果可逆，就删掉这一对。否则留下他（把 $x_i (b_i \wedge c_i)$ 再加回去）。
但是这样的话每次判断是否可逆的复杂度，即 Gauss 消元，是 $O(n^3)$，这肯定不行。怎么办呢？
如果是一般的情况，其实有 Sherman-Morrison-Woodbury 定理，可以做到每次 $O(n^2)$，总复杂度 $O(mn^2)$。不过本题中这个复杂度还是不行。
5. 最后的算法
我们发现其实在本题我们有一个很好的性质，那就是 $b_i$ 和 $c_i$ 最多只有最多两个分量非 $0$。所以 $x_i (b_i \wedge c_i)$ 只影响一个很小的子矩阵 $Y_{S, S}$。
因此我们可以使用 Harvey 定理——等一下，$Y$ 本来也不一定可逆欸。不过没关系，
回归正题，由于 $x_i (b_i \wedge c_i)$ 只影响一个很小的子矩阵，可以快速用 Harvey 定理判断更新后是否仍然满秩（可逆），如果可逆还可以快速求出更新后的逆。
不过这样，复杂度也还是 $O(mn^2)$。怎么进一步优化呢？
我们可以采用分治算法：如果我们有一堆边，他们对应的 $S$ （也就是他们连接的顶点对应的那些维）并起来是 $T$，而且这个 $T$ 不是很大的话，
我们就可以对他们统一更新：我们发现 Harvey 定理里判断去掉这些边之后是否合法只需要用到 $N_{T,T}$。所以我们可以只更新 $N_{T,T}$，不更新 $N_{*,*}$。
待到这些边都更新完之后，我们再去重新更新全部的 $N$。
如果把这个思路用到极致，就可以写出一个分治算法（伪代码）：
SOLVE(R, C):
记 S = R 并 C
// 执行时，要求 N_(S, S) 是正确的
// 执行后，已经尝试删除了所有一端点属于 R 而另一端点属于 C 的边
IF |R| > 2 or |C| > 2:
将 R 划分为 R1, R2
将 C 划分为 C1, C2
FOR (i, j) in { (1, 1), (1, 2), (2, 1), (2, 2) }:
SOLVE(Ri, Cj)
使用 Harvery 定理更新 N_(S, S)
ELSE IF R = { 2u-1, 2u }, C = { 2v-1, 2v } 且存在左端点为 u 右端点为 v 的边:
使用 Harvery 定理判断是否可以删掉这条边
如果是，删掉
// 不需要更新 N，因为递归上去之后会更新注意伪代码里默认了结点 $u$ 对应的行列还是 $2u-1, 2u$。事实上由于我们删了一些行列，可能未必如此。所以分治需要注意不要把一个点对应的两个列分开。
笔者的实现里是直接对结点分治，而非对矩阵行列分治。
这样的话，复杂度为 $T(n) = 4T(\frac{n}{2}) + O(n^w)$。根据主定理可以计算出 $T(n) = O(n^w)$。如果 $w = 2$，那么 $T(n) = O(n^2 \log n)$。
代码
由于这个矩阵是 $2n \times 2n$（如果大部分点都没有颜色并且答案几乎是满的），所以常数有点大。
笔者和出题人码风不通，没看懂他的实现。笔者本人的实现最后使用了 $O(n^3)$ 算法，不过模数取用了 $10^8 + 7$，这样可以在矩阵乘法时全部乘完加完最后取模，降低了很多常数。
#include <algorithm>
#include <cstdio>
#include <cstdlib>
#include <cstring>
typedef long long LL;
const int mod = 100000007;
const int N = 303;
const int NN = 606;
const int M = N * N / 2;
LL pow_mod(LL a, LL b) {
LL ans = 1;
for (; b; b >>= 1, a = a * a % mod)
if (b & 1) ans = ans * a % mod;
return ans;
}
inline void mul(int m, LL *A, LL x) {
for (int i = 0; i < m; ++i) A[i] = A[i] * x % mod;
}
inline void add(int m, LL *A, const LL *B, LL x) {
for (int i = 0; i < m; ++i) A[i] = (A[i] + x * B[i]) % mod;
}
inline void swap(int m, LL *A, LL *B) {
for (int i = 0; i < m; ++i) std::swap(A[i], B[i]);
}
void rank(int n, int m, LL A[NN][NN], int &r, int *S) {
// A[S][*] has rank r.
static int id[NN];
for (int i = 0; i < n; ++i) id[i] = i;
for (int i = 0, j = 0; j < m; ++j) {
int k = i;
while (k < n && A[k][j] == 0) ++k;
if (k == n) continue;
if (i != k) {
std::swap(id[i], id[k]);
swap(m, A[i], A[k]);
}
S[r++] = id[i];
LL t = pow_mod(A[i][j], mod - 2);
for (int k = i + 1; k < n; ++k)
add(m, A[k], A[i], -t * A[k][j] % mod);
++i;
}
}
bool inv(int n, const LL A_[NN][NN], LL B[NN][NN]) {
static LL A[NN][NN];
// return true iff A is invertible.
for (int i = 0; i < n; ++i)
for (int j = 0; j < n; ++j)
A[i][j] = A_[i][j], B[i][j] = i == j;
for (int i = 0; i < n; ++i) {
int k = i;
while (k < n && A[k][i] == 0) ++k;
if (k == n) return false;
if (i != k) {
swap(n, A[i], A[k]);
swap(n, B[i], B[k]);
}
LL t = pow_mod(A[i][i], mod - 2);
mul(n, A[i], t); mul(n, B[i], t);
for (int k = 0; k < n; ++k) if (k != i) {
add(n, B[k], B[i], -A[k][i] % mod);
add(n, A[k], A[i], -A[k][i] % mod);
}
}
return true;
}
void mul(int n, int m, int r, const LL A[NN][NN], const LL B[NN][NN], LL C[NN][NN]) {
static LL tmp[NN][NN];
for (int i = 0; i < n; ++i)
for (int j = 0; j < r; ++j) tmp[i][j] = 0;
for (int i = 0; i < n; ++i)
for (int j = 0; j < m; ++j)
for (int k = 0; k < r; ++k)
tmp[i][k] = tmp[i][k] + A[i][j] * B[j][k];
for (int i = 0; i < n; ++i)
for (int j = 0; j < r; ++j) C[i][j] = tmp[i][j] % mod;
}
#define take(A, X, Y, FORMULA) \
for (int i = 0; i < X; ++i) \
for (int j = 0; j < Y; ++j) \
A[i][j] = FORMULA
bool modify(const LL A[NN][NN], LL B[NN][NN], const LL A1[NN][NN],
int s, const int *S, int t, const int *T) {
// assuming B = A^(-1) and A[i][j] != A1 only when i, j in S
// if A1 is singular, return false
// else:
// modify B such that B_(T, T) = A1^(-1)_(T, T)
// only use A[S][S], A1[S][S], B[T][T].
static LL D[NN][NN], C[NN][NN], P[NN][NN];
take(D, s, s, A1[S[i]][S[j]] - A[S[i]][S[j]]);
take(P, s, s, B[S[i]][S[j]]);
mul(s, s, s, D, P, P);
for (int i = 0; i < s; ++i) (++P[i][i]) %= mod;
if (!inv(s, P, C)) return false;
mul(s, s, s, C, D, C);
take(P, t, s, B[T[i]][S[j]]);
mul(t, s, s, P, C, C);
take(P, s, t, B[S[i]][T[j]]);
mul(t, s, t, C, P, C);
for (int i = 0; i < t; ++i)
for (int j = 0; j < t; ++j) {
int ti = T[i], tj = T[j];
B[ti][tj] = (B[ti][tj] - C[i][j]) % mod;
}
return true;
}
int C[N], L[N];
LL X[N][N];
bool link[N][N];
int n, nn;
int id[NN];
int t, T[NN];
void Add(LL A[NN][NN], int u, int v, int i, int j, LL p) {
int x = L[u] + i, y = L[v] + j;
if (C[u] && i) x = L[u], p = p * C[u] % mod;
if (C[v] && j) y = L[v], p = p * C[v] % mod;
if (id[x] >= 0 && id[y] >= 0) A[id[x]][id[y]] = (A[id[x]][id[y]] + p) % mod;
}
inline void Add(LL A[NN][NN], int u, int v, LL p) {
// b = e(u0) - e(v0)
// c = e(u1) - e(v1)
// A += p (b wedge c)
Add(A, u, u, 0, 1, p);
Add(A, u, v, 0, 1, -p);
Add(A, v, u, 0, 1, -p);
Add(A, v, v, 0, 1, p);
Add(A, u, u, 1, 0, -p);
Add(A, u, v, 1, 0, p);
Add(A, v, u, 1, 0, p);
Add(A, v, v, 1, 0, -p);
}
LL Y[NN][NN], Z[NN][NN];
int ss[10], SS[10][NN];
LL AA[10][NN][NN], BB[10][NN][NN];
bool solve(int l1, int r1, int l2, int r2, int dep) {
#define copyS(U, V) \
for (int i = 0; i < s; ++i) \
for (int j = 0; j < s; ++j) \
U[S[i]][S[j]] = V[S[i]][S[j]]
if (r1 < l1 || r2 < l2) return false;
int &s = ss[dep] = 0;
int *S = SS[dep];
LL (*A)[NN] = AA[dep], (*B)[NN] = BB[dep];
for (int i = L[l1]; i < L[r1 + 1]; ++i) if (id[i] >= 0) S[s++] = id[i];
for (int i = L[l2]; i < L[r2 + 1]; ++i) if (id[i] >= 0) S[s++] = id[i];
std::sort(S, S + s); s = std::unique(S, S + s) - S;
copyS(A, Y);
copyS(B, Z);
if (r1 == l1 && r2 == l2) {
if (!link[l1][l2]) return false;
Add(A, l1, l2, -X[l1][l2]);
if (!modify(Y, B, A, s, S, 0, NULL)) return false;
link[l1][l2] = false;
Add(Y, l1, l2, -X[l1][l2]);
return true;
}
int m1 = (l1 + r1) / 2, m2 = (l2 + r2) / 2;
bool changed = false;
for (int a = 0; a < 2; ++a)
for (int b = 0; b < 2; ++b) {
int u1 = l1, v1 = m1, u2 = l2, v2 = m2;
if (a) u1 = m1 + 1, v1 = r1;
if (b) u2 = m2 + 1, v2 = r2;
if (solve(u1, v1, u2, v2, dep + 1)) {
changed = true;
// fprintf(stderr, "%d %d\n", ss[dep + 1], s);
if (ss[dep + 1] != s) {
if (false) {
inv(s, Y, Z);
} else {
copyS(Z, B);
modify(A, Z, Y, ss[dep + 1], SS[dep + 1], s, S);
}
}
copyS(A, Y);
copyS(B, Z);
}
}
return changed;
}
bool vis[N];
int path[N][N], len[N], ans;
bool dfs(int x, bool rt) {
vis[x] = true;
if (C[x] && !rt) {
++ans;
path[ans][len[ans]++] = x;
return true;
}
for (int y = 1; y <= n; ++y)
if ((link[x][y] || link[y][x]) && !vis[y]) {
int t = dfs(y, false);
if (t) {
path[ans][len[ans]++] = x;
return true;
}
}
return false;
}
int main() {
int m;
scanf("%d%d", &n, &m);
for (int i = 1; i <= n; ++i) {
scanf("%d", &C[i]);
L[i] = nn++;
if (!C[i]) nn++;
}
L[n + 1] = nn;
for (int i = 0; i < nn; ++i) id[i] = i;
for (int i = 0; i < m; ++i) {
int u, v;
scanf("%d%d", &u, &v);
if (u == v || link[u][v] || link[v][u]) continue;
link[u][v] = true;
X[u][v] = rand() % mod; // TODO: ERROR
Add(Y, u, v, X[u][v]);
}
rank(nn, nn, Y, t, T);
memset(id, -1, sizeof id);
for (int i = 0; i < t; ++i) id[T[i]] = i;
memset(Y, 0, sizeof Y);
for (int u = 1; u <= n; ++u)
for (int v = 1; v <= n; ++v)
if (link[u][v]) Add(Y, u, v, X[u][v]);
inv(nn, Y, Z);
solve(1, n, 1, n, 0);
for (int i = 1; i <= n; ++i)
if (!vis[i] && C[i]) dfs(i, true);
printf("%d\n", ans);
for (int i = 1; i <= ans; ++i) {
printf("%d", len[i]);
for (int j = 0; j < len[i]; ++j)
printf(" %d", path[i][j]);
puts("");
}
return 0;
}
鄂公网安备 42010202000505 号